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Algebraic algorithms for multiloop calculations 
The first 15 years. What's next? 

Fyodor V. Tkachov 
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The ideas behind the concept of algebraic ("integration-by-parts") algorithms for multiloop calculations are reviewed. For any 
topology and mass pattern, there exists a finite iterative algebraic procedure which transforms the corresponding Feynman- 
parameterized integrands into a form that is optimal for numerical integration, with all the poles in D-A exphcitly extracted. 



The Concept 

A victory lias a tiundred fathers 

The concept of algebraic algorithms for multiloop in- 
tegrals [l]-[4] has brought about an industry of analyti- 
cal multiloop calculations'", so it would not be inappro- 
priate to reminisce upon the anguish of its birth [1], 
given a non-negligible amount of research funds being 
misallocated through misinformation''. 

Back in 1979, when perturbative QCD, IBM, and the 
USSR were in full bloom, the first project which I took 
part in was successfully completed. Our analytical cal- 
culation was receiving a considerable attention, and my 
contribution to it was such that there was no doubt in my 
mind that it would be my Ph.D. thesis. However, the 
senior member of our team of three intimated to me that 
the result would go instead into another thesis that would 
not be able to get composed "otherwise" (which was true 
beyond any shade of doubt). The deal (in which I had no 
negotiation powers) was accompanied by a swift promo- 
tion of the senior member.'^ A concrete reward for me 
was a deterioration of my eyesight by two diopters as a 
result of my innocent calculational enthusiasm. Moreo- 
ver, due to the rigidity of Soviet system and my then in- 
significance, I was stuck in that vulnerable position for 
what seemed to be forever, although I was too light- 
hearted and inexperienced at the time to foresee all the 
consequences (which may have been just as well). 

But research was fun and the next level of complexity 



A number of papers in the proceedings of all AIHENP work- 
shops discuss calculations that make an essential use of identi- 
ties obtained via differentiations in momentum space. 

often undocumented; but cf. an unsuspecting citation in [5]. 

Of course, Germany (say) is rich enough to afford some over- 
head, especially if one takes into account the grant-fetching 
value of the Russian multiloop expertise. However, "tD()OSO is 
partner witl) a tl)ief l)atetl) l)is own souf: l)e l)earetl) curs* 
ing, an§ Besorayetfj it not." [6] 

'' Human nature is, of course, much the same everywhere in the 
world. But one ought to remember about behavioral and social 

renormalization factors — in the present case, roughly, 
Mr. James Carter over tovarishch Leonid Ilyich [Brezhnev]. 



(3-loop finite parts of massless self-energy diagrams) 
seemed within grasp. The appetites whetted, it was sug- 
gested that "we" tackle the many multiple series that 
emerged there within the technique that I had shown to 
be successful with 2-loop finite parts. However, I balked: 
if one were to appropriate the fruits of my meditations, 
one could at least leave me alone in the joy of the proc- 
ess. (That tacit arrangement was adhered to by the other 
side with much gusto all the way through the theory of 
asymptotic expansions, which is another story [7].) 

Being thus vigorously motivated, I conceived of a 
method that would, ideally, meet the following criteria; 

► My method would reduce human intervention to an 
absolute minimum. First, because my eyesight was at 
stake. Second, although I remember being able (then) to 
expand F-functions in 2-loops pretty reliably, I was less 
sure about 3 loops. 

► That meant computer algebra and Schoonschip [8], 
so that my ideal method would be basically simple to 
make a program implementation feasible. 

► Third, the usual approach was (and largely remains) 
to do all basic integrals "analytically" (i.e. by hand in the 
usual human fashion) and leave to computer only traces, 
scalar products, and substitutions of answers for the basic 
integrals. Since there were so many different integrals in 
3 loops, I decided that such an approach was unaccept- 
able: My method would have to do something about all 
those integrals too. 

► Lastly, my ideal method — my baby — had to be 
beautiful. That was the only kind of consolation I could 
expect to obtain with any degree of certainty. 

(NB The list is by no means obsolete after 15 years.) 

Where did I have to look for such a method? When- 
ever one has to start from scratch (which one has to do 
surprisingly often [7], [9]), a little philosophizing at the 
level of first principles is appropriate (because nothing 
works like a sound first principle implemented with due 
simplicity and relentless logic); 

One general consideration was not unknown (see e.g. 
[10]) but rarely heeded, namely, that best algorithms for 
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computer algebra applications should avoid imitating 
human ways of doing formulas. A rule of thumb is, one 
has to trade the true complexity of the amazing pattern 
matching that human beings are capable of, for a dull 
hugeness of more or less uniform data structures (e.g. 
polynomials) that computers are so good at. A difficulty 
here is that one has to understand quite well what's go- 
ing on inside computer algebra systems. Another diffi- 
culty is in imagining an "analytical" algorithm that one 
may not even be able to perform oneself. 

Furthermore, "calculation" is nothing but a "trans- 
formation". In order to "transform" an object one has to 
understand its "structure". What is "structure" of 3-loop 
massless self-energy integrals? (Exercise Answer this 
before reading on.) First of all, it is the set of identities 
expressing momentum conservation at vertices. That is 
already used in simplifications of scalar products in nu- 
merators. Next, "masslessness" means that Euler's iden- 
tity Qdf{Q) = Xf{Q) holds for massless self-energies 

with Q the external momentum. (At that point I already 
knew better to be bothered by comments about how it was 
silly to attempt to do integrals by differentiations.) Doing 
differentiation on the integrand and simplifying scalar 
products one obtains an identity connecting several dif- 
ferent integrals. The identity depends on how Q is passed 
through the propagators, so I started writing out all pos- 



sible variants and trying to combine them in all sorts of 
ways. It would have been easier to deal with 2-loop dia- 
grams, but I did not know the answer and chose the 3- 
loop ladder as my playing ground, so it took several 
months to stumble upon the celebrated triangle relation 
[1]. That took care of planar diagrams. A little more ef- 
fort yielded me non-planars in terms of one basic non- 
planar diagram (as well as a rebuke for having failed to 
explain sooner the entire significance of differential 
identities). I had a complete algorithm but could not re- 
sist the temptation of wasting some more time discover- 
ing (once a week or so) ways to reduce the one remaining 
scalar non-planar integral to planar ones before the algo- 
rithm was submitted for publication. 

But that was not the whole story. When actual pro- 
gramming of the algorithm based on pure recursions [1], 
[2] was attempted, it proved a failure (the huge program 
simply would not work as a whole). The project remained 
at a standstill until an explicit solution of the triangle 
recursion was found [3]. It resulted in a great speedup in 
both debugging and performance of the program which 
finally began to turn diagrams into numbers in a way that 
caused me to baptize it "Mincer". I emphasize a central 
role of that solution in the actual functioning of 
Mincer [4] because misleadingly few references have 
been made to the fact in the literature. 



What's Next? 

...For the duration of 40 hours he was expounding to 
Messrs. Gentlemen a great project, the execution of which 
subsequently brought a great glory to England and showed 
how far the human mind can reach out sometimes ! 
The Drilling Through The Moon With A Colossal Auger 
— such was the subject of Mr. Lund's speech !... [11] 

The range of applications of algebraic algorithms of 
the conventional type (differential identities in momen- 
tum space [1]) continues to expand (cf. the recent appli- 
cations to the electron magnetic moment [12] and to the 
effective heavy quark theories [13]). What I'd like to 
discuss now is the case of integrals with several external 
parameters (s, t, masses etc.), which is vital in modern 
high energy physics applications. 

First of all, with several parameters in the problem, 
the variety of analytical details is so enormous that one is 
driven by desperation and logic to ignore them all and to 
consider calculation of an arbitrary multiloop integral — 
however preposterous such an idea may seem! 

Second, one has to be clear about the meaning of 
"calculation". The enormity of the problem forces one to 
interpret it in a strictly pragmatic fashion as obtaining 
a family of "satisfactorily fast" algorithms to produce 
required curves. Such algorithms need not be based on 
"analytical" formulas. In fact, with many parameters in 
the problem, numerics at the final stage can hardly ever 
be avoided. 



Lastly, the true purpose of algebraic identities is not so 
much to "calculate" the integrals as to make them 
simpler OR reduce their number.'' 

NB The meaning of "simpler" is, again, strictly con- 
text-dependent: in the case of numerical integration it 
need not mean "fewer terms" in the integrand but e.g. 
"more continuous derivatives". For instance, the new 
multidimensional integration package Milx [15] 
provides an option of using higher-order quadratures 
along with Monte Carlo, thus greatly improving per- 
formance when integrands are sufficiently smooth. 

So, let x = ixi,...,Xj^) be the vector of Feynman pa- 
rameters; S is the integration region, a simplex described 
by Xj>0;l,jXj<l; j = l,...,K {K = dimSy, Q(x) 

and Vi (x) are arbitrary polynomials of x with symbolic 
coefficients; jj.- are symbolic (complex, transcendental) 
parameters. Consider integrals of the form 



Contrast this with the "algebraic calculability" of [14] which 

attempts to normalize the spectacular success achieved in the 
case of massless self-energies which was merely a matter of 
luck. A situation much more like what I envisaged in the be- 
ginning occured in the case of anomalous magnetic moment 
where everything reduces to 18 independent integrals [12]. 
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jdxQix)UVi^'(x). (1) 

s ' 

Dimensionally/analytically regulated diagrams are ob- 
tained by appropriate choices of Qix) , V) (x) and . 

Hi are equal to a fixed integer (denoted n,- ; it can be 
positive, negative or zero) plus a symbolic part (denoted 
fj ): £j are regarded as infinitesimal because an expan- 
sion in £■ has to be performed ultimately. The point is 
that, whereas the integral (1) may be well-defined as a 
meromorphic function of /i,- , the expansion in e- may 
not be performed directly in the integrand because some 
rii may be too large and negative. 

The structure I propose to exploit can be described as 
"polynomials to symbolic powers" and is fully specified 
by (1): no further information about V, (jc) will be used. 

The key result I would like to present is as follows: 

Theorem A ("algorithm of algorithms of algebraic 
multiloop calculations"). Any integral (1) can be re- 
duced — via a finite number of purely algebraic steps — 
to the following form with any given integer oa > .• 

£ ^\^' a{x)X\Vr''ix), (2) 

where the sum contains a finite number of terms; S' can 
be either S or any of its boundary simplices with vari- 
ous dimS'<S (dim5"=0 is allowed); x parameterize 
5" in a natural fashion (e.g. some sets x are obtained 
by nullifying certain components of x); Q'(x') are 

polynomials of x ; B' are independent of x ; B' and all 

coefficients of Q'{x) are polynomials of all Hi and of 
the coefficients of Q{x) and y,(x); Vjix ) are restric- 
tions of Viix) to the corresponding boundary simplices 
S' ; all numeric coefficients involved are integer. 

The proof uses the following result: 

Theorem B ("generalized Bernstein functional 
equation"). For any finite set of polynomials Vi{x) 

there exists an identity of the following form: 

B-ifp {x,d)Y[ v^i^^ ix) = n v;.^' (x) , (3) 

where '2 (x,d) is a polynomial of x and di = d I dxi ; B 
and all coefficients of '£ are polynomials of Hi ^nd of 
the coefiicients of V- (x) . 

Proof of Theorem B Using methods of abstract 

algebra ("finitely generated ideals" of "graded left mod- 
ules", etc.), one can establish (see e.g. [16]) existence of 
identities (3) for the case of exactly one (arbitrary) poly- 
nomial — a result known as the Bernstein functional 
equation: 



b-^!P(x,d)Vf+\x) = Vf(x), (4) 

where both b and !P depend polynomially on fi and on 
the coefficients of V. Take y (x) = 0, Vi (x) , perform the 
differentiations in (4), cancel the symbolic powers and 
equate coefficients of monomials of x on both sides; one 
obtains a linear system of equations for the unknown 
coefficients of (P . Bernstein's result is equivalent to ex- 
istence of a non-zero solution of the system, i.e. some 
determinant is non-zero. Perform a similar reduction for 
(3), and compare the resulting linear system with that 
obtained for (4). The former differs from the latter by 
replacements of different entries of fi by different Hi > so 
that the determinant mentioned above cannot be zero. 
Therefore, the system has a non-zero solution, which 
yields the identity (3) with the coefficients of fP that are 
rational functions of both Hi the coefficients of 
Vi (x) ■ Taking out all the denominators and combining 
them into B yields (3). □ 

Proof of Theorem A One iteratively performs re- 
placements of the integrand with the corresponding l.h.s. 
of (3) and performs integrations by parts to get rid of the 
differentiations. One does this (with appropriate rp and 
B) also for all integrals corresponding to the boundary 
terms resulting from integration by parts. After a finite 
number of iterations the integer parts of all complex 
powers become > to . □ 

In dimensionally regulated diagrams, all are linear 
combinations of one symbolic parameter s. This may 
cause problems due to nullification of some B. That 
would mean a breakdown of dimensional regularization 
so that an additional analytic regularization would have 
to be used to make all jUi independent. But since the 
identities (3) have to be found for arbitrary anyway (to 
handle integrals with different integer parts of /<;•), there 
is no practical problem here. ■ 

Discussion 

• The operators fP for the same topology but different 

mass patterns need not be connected in any simple way. 

• There is an infinite family of non-trivially different 
operators fP satisfying (3) for the same V) (x) . 

• The only way, known to me at the time of this writ- 
ing, of finding the polynomials T is via a direct study of 
the linear system of equations for the unknown coeffi- 
cients of fP (cf. the proof of Theorem B). This is cumber- 
some (the form of the system depends on the degree of rp 
w.r.t. derivatives), but one has to find such an identity 
only once for each topology/mass pattern, and only poly- 
nomial algebra with rational coefficients is involved. 

• Multiplying both sides of (3) from the left by some of 
V, (x) 's and absorbing the latter into rp on the l.h.s., one 
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proves existence of identities with only a subset of expo- 
nents raised by +1, with others intact. In general, the 
family of operators of this latter form contains ones that 
cannot be obtained from (3) as just described. These ad- 
ditional identities may result in more economical algo- 
rithms (iteratively raising each symbolic power to exactly 
the minimal value required for achieving the desired 
smoothness of the integrand). 

• Similarly, one can directly seek an identity wdth -1-1' s 
on the l.h.s. replaced with + N^'s with required Ni > 
(or even with some A^,- < if the integer part of the i'-th 

original power was already positive), and with Q on the 
r.h.s. This is similar to explicit solution of recurrent re- 
lations in the conventional algorithms [3], [12]. 

• More iterations make integrands smoother, but the 
original integrals did have singularities in external pa- 
rameters! The paradox is resolved by noticing that the 
denominators B have zeros. Those zeros, therefore, must 
correspond to both threshold singularities (hence B must 
be connected to determinants of the system of Landau 
equations) and to the poles in D-4 (in the case of di- 
mensional regularization). 

• It follows that the algorithm automatically extracts 
the poles in D-4 : just expand integrands in the end in 
D-4 ; complex powers give rise to integer > powers 

of Viix) and logl/ (x) ; all the poles come from B's. 
With well-defined integrals for the coefficients of the 
poles, their cancellations may be checked numerically. 

• Due to zeros of B's, the numerical advantage of 
smoother integrands is lost near thresholds (where nu- 
merical convergence is poor anyway) due to cancellations 
in the sum (2). However, the regions near thresholds may 
be systematically treated via asymptotic expansions ob- 
tained by the method of asymptotic operation [7]. 

• The identities of the original "integration-by-parts" 
algorithm [l]-[2] can be regarded as special cases of (3). 
In that sense, Eq.(3) represents the most general basic 
type of algebraic identities for multiloop integrals, which 
justifies the name of Theorem A. Concerning the identi- 
ties connecting integrals with different values of the 
space-time dimension D (cf. [18], [20], [19]), one notices 
that D enters exponents of different factors with dijferent 
signs, whereas in (3) all exponents change in the same 
direction. So one cannot naturally interpret the new 
identities as recurrences w.r.t. D (beyond one loop). 

• As a simple application, consider dimensionally 
regulated one-loop integrals which contain just one com- 
plex exponent, Jj dx Q(x)y~"~^(jc) , where n is a non- 
negative integer, e is the complex regularization pa- 
rameter, and V{x) is a quadratic polynomial of x. 

Without loss of generality V{x) = x^V x+2R^ x + Z 

where V is a KxK matrix and 7? is a AT-vector. Then 



1- 



{x+A)d, 



2(M + 1) 



V^'+\x) = V^'ix), 



(5) 



where A : 



^Rj and A = R^V-^. (Special 



cases of this relation in integral form have been known 
for some time; cf. [18].) (Exercise Write this out for 
your favorite one-loop integral.) (NB If V is singular 

then one should determine its null space and treat trans- 
verse and longitudinal coordinates separately.) Usually, 
one first reduces non-scalar one-loop integrals to scalar 
ones and then uses formulas for the latter in terms of 
dilogarithms etc. — a typical example of a procedure 
motivated by how human beings perform calculations. 
The computer-oriented approach based on (5) is to apply 
it iteratively to the integrand (with all the numerators) 
until required regularity is achieved; then one can take 
advantage of the numerical integration schemes that em- 
ploy higher-order quadratures such as Milk [15].® 

• For 2-loop diagrams, the situation is as follows. For 
the simplest sunset topology, if !P is second order w.r.t. 
d's (a minimal possible number of derivatives that 
yields a non-trivial identity (3)) then has 5-th powers 
of x's. (NB The maximal powers of d's and x's are anti- 
correlated.) This is quite cumbersome, but it is a dull 
complexity of a polynomial. 

• Being a monotonous repetition of the same cumber- 
some algebraic identity, the algorithm is as unfitted for 
hand calculations as it is perfect for programming. 

Conclusions 

It took about 10 years for Mincer [4], an implementa- 
tion of the first algebraic multiloop algorithm [l]-[3], to 
deliver reliable results. The amount of algebra involved 
in the construction and iterative application of (3) is even 
greater, and may prove huge enough to prevent the new 
"algorithm of algorithms" from being actually used be- 
yond one loop to any greater extent than the SONG OF 
Songs is sung. 

But consider this: If an identity (3) is explicitly con- 
structed for a given topology/mass pattern, then it simul- 
taneously does the following things — each of which 
would normally be an achievement per se: 

^ It reduces integrands to a form that is optimal for 

numerical integration. 

^ It allows straightforward program implementation. 

^ It explicitly extracts all (!) the poles in D-4 . 

Taking into account the unrestricted generality of the 
new approach (it is applicable — at least in theory — 
to any multiloop integral, including ones with, say, non- 
relativistic propagators etc. etc.) and the pace of progress 



' A specialized version of Milx is to be made available. 
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of computer industry, it looks like it might possess a po- 
tential to relieve us from some grueling chores in future 
and to give a boost to the industry of loop calculations 
(which needs it [21]; remember also the funds that went 
into the lattice QCD). The operators rP could be found by 
solving a system of linear algebraic equations for their 
coefficients — a highly cumbersome procedure given all 
the symbolic parameters involved {D, masses etc.) but 
requiring nothing that would go much beyond the stan- 
dard arsenal of polynomial computer algebra. So, anyone 
equipped with a sufficiently fast and reasonably flexible 
computer algebra system may already try one's hand at 



constructing an identity (3) for one's favorite two-loop 
topology. 
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